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ABSTRACT 


A  two-dimensional  plasma  actuator  analysis  code  has  been  developed.  A  time-accurate  Navier-Stokes 
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in  quiescent  air  over  a  flat  plate  using  only  the  plasma  actuator  model.  The  resulting  fluid  motion  was  similar  to 
the  induced  motion  observed  in  plasma  actuator  experiments. 
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INTRODUCTION 


Aerodynamic  flow  control  by  means  of  plasma  actuators  is  an  active  area  of  research  in  the  Air  Force 
and  in  academia.  Experiments  have  shown  that  these  devices  can  be  used  to  induce  flow  movement  in  a 
stationary  air  mass,  and  to  reattach  separated  flow  on  airfoils  at  high  angles  of  attack.  Plasma  actuators  offer 
potential  advantages  in  weight  reduction,  cost,  and  manufacturability  over  traditional  flow  control  devices  such 
as  slats  and  flaps1. 

Plasma  actuators  consist  of  two  electrodes  separated  by  a  dielectric  material.  The  electrodes  are  offset, 
with  one  electrode  exposed  to  the  air  and  the  other  embedded  in  the  dielectric  material  as  shown  in  Figure  1. 
When  a  high  voltage  alternating  current  is  applied  to  the  electrodes,  the  air  above  the  insulated  electrode 
ionizes,  beginning  at  the  edge  of  the  exposed  electrode  and  spreading  across  the  insulated  electrode.  The  extent 
to  which  the  plasma  field  covers  the  insulated  electrode  increases  and  decreases  with  the  rising  and  falling 
voltage.  The  ionized  region  produces  a  “body  force”,  which  induces  motion  in  the  surrounding  air2. 


Figure  1:  Plasma  Actuator  Schematic3 

In  FY  2004-2005,  AFOSR  sponsored  AFRL/MNAC  to  perform  research  into  the  effectiveness  of 
plasma  actuators  as  aerodynamic  control  devices.  The  ultimate  objective  of  this  research  was  to  model  the 
effects  of  plasma  actuators  on  aerodynamic  flow  fields  by  incorporating  a  first-principles  plasma  generation 
model  into  a  Navier-Stokes  flow  solver.  Both  the  plasma  code  and  the  flow  solver  code  were  under 
development  in  FY  2004-2005.  While  the  flow  solver  had  been  successfully  demonstrated  on  several  standard 
validation  cases,  development  of  the  plasma  code  had  not  progressed  to  the  point  that  it  could  be  incorporated 
into  the  flow  solver.  Because  of  this,  AFRL/MNAC  initiated  a  collaborative  agreement  with  AFOSR-sponsored 
researchers  at  the  University  of  Notre  Dame  to  incorporate  a  plasma  actuator  model  developed  there  into  the 
AFRL/MNAC  flow  solver. 

The  Notre  Dame  model  may  best  be  described  as  phenomenological.  While  this  model  does  not 
simulate  the  creation  and  motion  of  ions  in  the  plasma  field,  it  does  a  credible  job  generating  force  fields  that 
induce  the  type  of  motion  observed  in  plasma  actuator  experiments.  Instead  of  simulating  the  plasma  field,  this 
model  represents  the  region  over  the  insulated  electrode  as  a  set  of  parallel  electronic  circuits,  such  as  the  one 
shown  in  Figure  2,  distributed  along  the  insulated  electrode.  The  circuit  capacitances  and  resistances  are  derived 
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from  both  experimentally  determined  and  inferred  electrical  properties  of  air.  The  individual  circuits  are 
activated  when  the  voltage  difference  between  the  exposed  and  insulated  electrodes  exceeds  a  critical  threshold. 
In  this  way  the  advance  and  retreat  of  the  plasma  field  is  simulated  as  the  applied  voltage  oscillates. 

—  Exposed  electrode 

—  Dielectric  surface 


Insulated  electrode 

Figure  2:  Plasma  Actuator  Circuit  Model2 

The  voltage  distribution  resulting  from  the  solution  of  the  time-varying  circuit  equations  is  used  on  the 
boundary  of  a  two-dimensional  Poisson  equation  solver  that  computes  the  electric  potential  in  the  vicinity  of  the 
plasma  actuator.  A  body  force  field  calculated  from  this  electric  potential  is  used  in  the  Navier-Stokes  flow 
solver  to  induce  motion  in  the  flow. 

The  remainder  of  this  report  will  describe  the  solution  algorithms  employed  in  the  two-dimensional 
Navier-Stokes  flow  solver  with  the  integrated  plasma  actuator  model.  Several  example  solutions  will  be 
presented  to  demonstrate  the  flow  solver’s  capability. 

TECHNICAL  APPROACH 

The  flow  solver  developed  in  this  effort  solves  the  Navier-Stokes  equations  of  fluid  motion  using 
simple  numerical  algorithms.  Although  solution  times  could  be  reduced  by  using  more  complicated  solution 
techniques,  simple  algorithms  were  chosen  in  order  to  reduce  the  effort  required  to  understand  the  code  well 
enough  to  perform  modifications. 

Since  the  purpose  of  this  research  effort  was  to  study  the  effects  of  plasma  actuators  on  air  flows, 
extraordinarily  high  solution  fidelity  was  not  a  requirement.  The  flow  solver  was  only  required  to  produce 
solutions  of  sufficient  fidelity  that  the  effects  of  the  plasma  actuator  on  the  flowfield  could  be  studied. 

In  this  finite-volume  flow  solver,  the  spatially  varying  terms  are  discretized  using  central  differencing4 
stabilized  with  flux-limiting  numerical  dissipation5.  Turbulent  terms  are  computed  using  the  so-called 
“detached  eddy  simulation”  version  of  the  Spalart-Allmaras  turbulence  model'1.  Temporal  integration  is 
performed  using  the  dual-time  stepping  algorithm7.  The  plasma  actuator  is  simulated  using  the  Notre  Dame 
model2'8. 
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Equations  of  Fluid  Motion 

The  two-dimensional  Navier-Stokes  equations  of  fluid  motion  in  conservation  law  form,  accounting 
for  body  forces,  are9 


5Q  dE  dF  _ 
dt  dx  dy 


5Q  t  a(E,-Ev)  (  8(Fj  -Fv) 
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The  shear  stresses  (t»*)  and  heat  fluxes  (q*)  in  Eq  (3)  are  given  by: 
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The  gas  properties  appearing  in  Eqs  (3)-(4)  are: 
p  Density 

u,v  Cartesian  velocity  components 
E  Total  energy 

fx,  fy  Cartesian  body  force  components  (per  unit  volume) 
p  Pressure 

H  Total  enthalpy  (E  +  p/p) 

T  Temperature 

|i  Coefficient  of  viscosity 

k  Coefficient  of  thermal  conductivity 

Two  equations  of  state  can  be  used  to  relate  the  fluid  properties: 


(4) 
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p  =  (y  - 1)  pE  -  -jp(u2  +  v2 )  or  p  =  pRT 


(5) 


Finite  Volume  Formulation 

The  Navier-Stokes  equations  are  solved  using  a  cell-centered  finite  volume  (FV)  formulation.  In  the 
FV  scheme,  the  domain  of  interest  is  subdivided  into  small  control  volumes,  or  cells.  A  portion  of  a  two- 
dimensional  computational  grid  is  shown  in  Figure  3.  The  details  of  this  diagram  will  be  explained  in  the 
following  discussion. 
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Figure  3:  A  Portion  of  a  Two-Dimensional  Computational  Grid 

The  fluid  properties  are  assumed  to  be  constant  throughout  a  particular  grid  cell  at  any  time  during  the 
solution  process.  Temporal  changes  in  the  cell  fluid  properties  are  determined  by  tracking  the  flux  of  the  fluid 
properties  across  the  cell  boundaries.  Eq  (1)  can  be  expressed  as 

^  +  V.(Ei  +  Fj)=S  (6) 

Treating  the  fluid  properties  within  a  cell  as  constant  and  integrating  Eq  (6)  over  the  cell  volume  V  yields 

vf+JJJv.(Er  +  FjK  =  vs  (7) 

Vol 

where  At*)  indicates  a  discrete  change  in  the  indicated  quantity.  Gauss’  theorem  can  be  used  to  convert  the 
volume  integral  of  Eq  (7)  into  a  surface  integral: 

V77  +  S(eJ  +  F^)  -"nAn-VS  (8) 

At  N=1  N 
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In  Eq  (8)  the  summation  is  over  the  bounding  faces  of  the  cell,  which  number  Nf.  AN  is  the  area  of  bounding 
face  N,  and  nN  =  nxi  +  nyj  is  an  outward  directed  unit  vector  normal  to  face  N  (see  Figure  1). 

The  inviscid  and  viscous  terms  can  be  treated  separately: 

v^  +  Ri(Q)-Rv(Q)- vs=  0  (9) 

where  the  inviscid  and  viscous  residuals  are 

Ri(Q)=^(Eii  +  F1j)N.]iNAN  (10) 

N=1 

Nf  /  x 

Rv(Q)=XlEvi  +  FvJ/N  '"nAN  (H) 

N=1 


Central  Differencing 

In  the  central  difference  scheme,  Eq  (10)  can  be  written 


Nf 

RI(Q)  =  XF(Q)NAf 


(12) 


where 


F(Q)n  —  (^N  '  ^n)Qn  + 


Pn 


0 

nx 

_  ny 
VN  •  n. 


(13) 


In  Eq  (13),  subscripted  flow  properties  (*)N  represent  the  average  value  of  the  indicated  flow  property  at  cell 
face  N.  For  example, 


VN=|(v0+Vn) 


(14) 


where  V0  is  the  velocity  in  the  cell  under  consideration  and  Vn  is  the  velocity  in  the  cell  adjacent  to  cell  0  at 
face  N  (see  Figure  1). 

The  viscous  fluxes  involve  partial  derivatives  of  the  Cartesian  velocity  components  and  the 
temperature.  The  partial  derivatives  for  a  fluid  property  in  cell  0  can  be  calculated  from10 

Nf 


ff1)  -f 

V  <3x  Jo  V0  ^ 

A-h  .2-2(4^ 


(15) 


dy 


0  V0  N=1 


Here,  (•)  represents  the  property  (e.g.  T)  to  be  differentiated. 

The  coefficient  of  viscosity  is  calculated  from  Sutherland’s  law 
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where 


Cj  =  1.458  xlO“ 


kg 


m  -  sV°K 


,C2  =110.4°K 


(17) 


The  coefficient  of  thermal  conductivity  is  calculated  from 

k_jL_ 

(y-l)Pr 


(18) 


where  y  is  the  ratio  of  specific  heats  (1.4  for  air),  R  is  the  gas  constant  (287  m2/s2  °K  for  air),  and  Pr  is  the 
laminar  Prandtl  number  (0.72  for  air). 

The  velocity  and  temperature  gradient  terms  are  calculated  for  each  cell  using  Eqs  (15).  Then  the 
viscous  residual  term  can  be  calculated  using  Eq  (1 1)  where,  for  example 


Artificial  Dissipation 

Flow  solutions  calculated  using  central  difference  spatial  discretization  schemes  sometimes  exhibit  a 
numerical  phenomenon  known  as  odd-even  decoupling,  in  which  alternating  grid  points  converge  to  different 
solutions.  In  order  to  stabilize  the  calculations  when  solving  the  Euler  or  Navier-Stokes  equations  using  a 
central  difference  scheme,  it  is  often  necessary  to  add  a  dissipative  term  to  the  FV  integral.  The  numerical 
dissipation  scheme  of  Yoon  and  Kwak5  is  used  here.  This  scheme  provides  stabilizing  background  dissipation 
to  the  central  difference  discretization,  and  employs  flux  limiters  to  smooth  out  numerical  dispersion  around 
sharp  gradients  such  as  shocks. 

At  each  cell  interface  a  dissipative  term  is  added  to  the  flux.  At  a  cell  interface  in  the  i-direction  of  a 
structured  grid,  for  example,  the  dissipative  term  is 


(20) 


In  Eq  (20),  the  (*)i+i/2  terms  are  evaluated  at  the  cell  faces  and  the  (#)j  terms  are  evaluated  at  the  cell  centers.  The 
coefficient  is: 


+  a 


'i+l/2 


(21) 


where  Kq  provides  a  threshold  dissipation  and  Kj  ensures  adequate  dissipation  in  the  vicinity  of  a  shock.  Ai+i/2  is 
the  cell  face  area  and  vi+1/2  is  a  “pressure  switch"  that  becomes  large  in  the  vicinity  of  a  shock: 


(22) 


vi+i/2  =max(vi_i,Vi,vi+1,vi+2) 

v  |Pi+i  -2Pi  +  Pi-i  | 
(Pi+i+2pi+pi_1) 

(|)  and  \|/  are  flux-limiting  functions: 

if  c  <  0 
if  0  <  ct  <  1 
if  1  <  a 

Y|/(a)  =  40  j 

where 


and,  for  example,  when  computing  dissipation  for  the  continuity  equation: 

ei .+y2  =  Pw  _  Pi 


(23) 

(24) 

(25) 

(26) 

(27) 


Explicit  Time  Integration 

Eq  (9)  can  be  rewritten  to  indicate  numerical  integration  over  a  discrete  time  interval: 

Qn+1=Qn  _  AtR(Q)n  (28) 

where  R(Q)  =  R,(Q)  -  R/Q)  -  VS  and  the  superscripts  (*)n  and  (*)n+1  indicate  flowfield  properties  evaluated  at 
the  current  solution  time  t  and  at  the  next  solution  time  t+At,  respectively. 

Eq  (28)  is  the  single-step  Euler  explicit  time  integration  from  time  level  n  to  time  level  n+1  and  has 
first  order  temporal  accuracy.  Multistage  explicit  integration  schemes  have  greater  temporal  accuracy  than  this 
single  stage  scheme.  For  example,  an  M-stage  Runge-Kutta  scheme  has  a  formal  accuracy  of  M  for  a  linear 
equation11.  In  the  M- stage  Runge-Kutta  scheme,  the  solution  is  advanced  over  the  time  interval  At  by 

q(°)  =  Qn 

Q(m)  =  Q(o)  — A.(m)^-R(Q(m-i))?  m  =  1,2,... ,M  (29) 

Qn+1  =  Q<m) 

Here,  the  parenthesized  superscripts  denote  the  stages  of  the  Runge-Kutta  scheme.  The  weighting  factors  X(m) 
are  calculated  from 

A-(m)  =  — ^ - ,  m  =  l,2,...,M  (30) 

M  +  l-m 
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The  allowable  time  step  for  stable  calculations  is  a  function  of  wave  propagation  speed  through  the 
grid  cells.  For  time-accurate  solutions,  the  time  step  At  is  the  global  minimum  of" 


At  < 


(Fx+Fy) 


(31) 


with 

fx={u|  +  c)Ax 

Fy={v|  +  c)Ay 


(32) 


where  V  is  a  cell  volume,  C  is  the  speed  of  sound  in  the  cell,  Ax  and  Ay  are  the  projected  areas  of  the  cell  in  the 
x-  and  y-directions.  Multistage  Runge-Kutta  integration  allows  the  time  restriction  to  be  relaxed,  so  the  time 
step  calculated  using  Eq  (31)  can  be  multiplied  by  a  factor  CFL  >  1. 


Accelerating  Convergence 

Convergence  to  a  steady  state  solution  can  be  accelerated  using  local  time-stepping  and  implicit 
residual  smoothing  during  the  Runge-Kutta  integration11. 

Local  time-stepping  dispenses  with  global  time  accuracy  by  advancing  the  solution  in  each  cell  using 
the  time  step  calculated  for  that  cell  using  Eq  (31).  This  drives  the  global  solution  error  rapidly  to  zero  since  the 
solution  errors  are  propagated  through  each  cell,  and  hence  through  and  out  of  the  solution  domain,  at  the 
maximum  possible  rate. 

Implicit  residual  smoothing  gives  the  explicit  time  integration  scheme  an  implicit  character,  allowing 
solution  advancement  at  time  steps  greater  than  indicated  by  the  stability  limit  of  the  explicit  scheme.  A 
Laplacian  smoothing  operator  is  applied: 

Rq  =  Ro  +  eV“Rn  (33) 


where 


Nf 

V2RN=£(Rn-R0) 


N=1 


Eq  (33)  is  solved  using  Jacobi  iteration: 


TJ  (°)  _  R  n 
^0  _  ^0 


rM  _  a_ 


Ro’+e^X 


Nf 

5"  (m-0 


N=1 


Nf 


m  =  1,2,... ,m 


MAX 


1  +  E^1 


N=1 


Ro  =  R0 


(34) 


(35) 
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Effective  values  of  the  weighting  factor  s  range  from  0.5  to  0.8.  In  Reference  11,  the  authors  reported  a 
doubling  of  the  allowable  time  step  for  Euler  solutions  on  unstructured  grids  by  using  two  Jacobi  iterations  per 
Runge-Kutta  stage  with  s  =  0.5. 


Dual  Time-Stepping 

Dual  time-stepping  (DTS)  is  an  implicit  time  integration  scheme7.  In  DTS,  the  spatial  terms  of  Eq  (6) 
are  discretized  at  time  n+1,  rather  than  at  time  n,  and  the  temporal  term  is  discretized  with  a  three-point 
backward  difference,  resulting  in  the  following  variation  of  Eq  (9): 

3Qn+1  -4Qn  +Qr 


V 


2At 


=  -R 


(q"+1  ) 


(36) 


This  equation  is  solved  by  introducing  Q  ,  an  approximation  to  Qn+  ,  and  a  pseudo-time  variable  t  . 
The  following  equation  results: 


V^  =  -R 
dt 


(Q-) 


-v 


3Q*  -  4Qn  +  Q 
2  At 


n-A 


=  -R 


(Q-) 


(37) 


Eq  (37)  can  be  solved  using  the  explicit  Runge-Kutta  scheme  of  Eq  (29)  with  local  time-stepping  in 
pseudo-time  (t  )  and  implicit  residual  smoothing.  Time  accuracy  is  ensured  by  the  inclusion  of  the  time 
derivative  term  in  R*(  Q*).  Eq  (37)  is  solved  iteratively  until  it  converges  on  a  solution  for  Q*: 

1.  Set  k  =  0 

2.  Set  Qk  =  Qn 

3.  Repeat  until  Qk+1  ~  Qk: 

a.  Solve  Eq  (37)  for  Qk+1  using  Eq  (29)  with  local  time-stepping  and  implicit 
residual  smoothing,  Q(0)  =  Qk,  and  Qk+1  =  Q(M) 

b.  Set  k  =  k  +  1 

4.  Set  Qn+1  =  Qk 


Turbulence  Model 

Inclusion  of  a  turbulence  model  was  not  originally  planned  by  the  MNAC  researchers  because 
turbulence  was  not  expected  to  play  a  significant  role  in  the  flows  to  be  investigated.  However,  as  the  code  was 
being  evaluated  against  standard  CFD  code  validation  cases,  it  became  apparent  that  unsteady,  separated 
flowfields  were  developing  in  laminar  solutions  of  test  cases  that  were  known  not  to  exhibit  such  behavior.  It 
was  decided  to  incorporate  a  turbulence  model  into  the  flow  solver  to  ensure  that  the  code  would  produce 
separated  flows  only  in  those  cases  where  it  could  reasonably  be  expected  to  occur. 

Turbulent  terms  are  computed  using  a  variation  of  the  Spalart-Allmaras  one-equation  turbulence 
model.  The  Spalart-Allmaras  equation,  written  without  transition  terms,  is7: 

1^  +  v(vv)=CblSv  +  i{v[(vL  +v)Vv]  +  Cb2(VvXVv)}-Cwlfw^j  (38) 
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Eq  (38)  is  solved  forV  ,  from  which  the  turbulent  viscosity  coefficient  can  be  calculated.  The  various  terms 
appearing  here  are  defined  in  References  6  and  7.  The  variation  from  the  standard  Spalart-Allmaras  model  lies 
in  the  d  term  which  appears  explicitly  in  Eq  (38)  as  well  as  implicitly  in  the  definitions  of  several  of  the  other 
terms  in  the  equation.  In  the  original  scheme,  d  is  the  distance  from  a  cell  center  point  to  the  nearest  solid  wall 
boundary  point.  Here,  d  is  the  smaller  of  this  distance  and  the  local  grid  spacing.  In  the  near  vicinity  of  a  solid 
wall,  Eq  (38)  behaves  as  a  turbulence  model,  but  away  from  the  wall,  it  assumes  the  nature  of  the  subgrid  scale 
closure  models  used  in  large  eddy  simulation  codes. 

The  solution  of  Eq  (38)  will  not  be  discussed  in  great  detail  here.  The  spatial  discretization  is  a  first 
order  upwind  scheme7,  and  the  temporal  solver  is  the  multistage  Runge-Kutta  scheme  of  Eq  (29)  with  local 
time-stepping.  Because  Eq  (38)  depends  on  the  current  flow  field,  it  is  not  solved  simultaneously  with  the 
Navier-Stokes  equations.  Rather,  Eq  (38)  is  advanced  one  local  time  step  for  each  stage  of  the  flowfield  Runge- 
Kutta  integration.  Using  this  approach,  Eq  (38)  can  be  driven  to  a  steady-state  solution  in  the  limit  of  a  steady- 
state  flowfield  solution. 

The  turbulent  viscosity  coefficient  p,  is  calculated  fromv  ,  and  is  added  to  the  laminar  viscosity 
coefficient  calculated  using  Eq  (16),  leading  to  the  following  redefinition  of  p: 

fi  =  fii+Pt  (39) 


The  coefficient  of  thermal  conductivity  for  turbulent  flow  is  calculated  by 

f-  --  \ 


k  = 


yR 


(r-i) 


Hi 

Pr 


JH_ 

Pr, 


(40) 


Prt  is  the  turbulent  Prandtl  number  (0.9  for  air). 


Plasma  Actuator  Model 

The  plasma  actuator  model  provided  by  the  University  of  Notre  Dame  uses  a  distribution  of  A  parallel 
electronic  circuits,  as  shown  in  Figure  3,  to  simulate  the  response  of  the  region  above  the  insulated  electrode  to 
an  applied  voltage8.  Each  subcircuit  represents  a  portion  of  the  insulated  electrode  length  and  has  finite  width 
and  length.  The  first  subcircuit  has  the  shortest  length  while  the  N,h  subcircuit  has  the  longest.  Each  subcircuit 
has  two  capacitors:  one  for  the  air  above  the  dielectric  and  one  for  the  dielectric  layer.  The  plasma  resistance  is 
modeled  with  different  values  for  the  forward  (positive  plasma  current)  and  backward  (negative  plasma  current) 
portions  of  the  AC  cycle.  Diodes  turn  the  resistance  subcircuits  on  in  the  presence  of  plasma  and  off  in  its 
absence. 
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Figure  3:  Parallel  Plasma  Actuator  Circuit  Model8 


The  air  capacitance  is  defined  as: 


Cm  = 


£o£aAn 

L 


(41) 


where  s0  is  the  permittivity  of  free  space  (8.854><10"12  F/m),  ea  is  the  dielectric  coefficient  of  air,  An  is  the  cross- 
sectional  area  of  the  air  capacitor,  and  ln  is  distance  from  the  edge  of  the  exposed  electrode  to  solution  point  n 
on  the  dielectric  surface. 

The  dielectric  capacitance  is: 


Cdn  = 


£0£d  Ad 

la 


(42) 


Here,  e(1  is  the  dielectric  coefficient  of  the  dielectric  material,  A(1  is  the  cross-sectional  area  of  the  dielectric 
capacitor,  and  ld  is  the  thickness  of  the  dielectric  material. 

The  subcircuit  resistance  is: 


R 


n 


PA 

A„ 


(43) 


where  pa  is  the  resistivity  of  the  air  and  An  and  ln  are  the  cross-sectional  area  and  length  of  the  subcircuit. 
Different  values  are  assigned  to  R„  depending  on  the  direction  of  current  flow;  Rn,  is  used  to  denote  the 
resistance  for  forward  plasma  current  flow,  while  Rnb  is  used  for  backward  current  flow. 

Given  a  time-dependent  applied  voltage,  the  voltage  on  the  dielectric  surface  in  subcircuit  n  is 
described  by 


dVn(t)  dVapp(t), 

(  c  1 

'-'an 

+  kn 

r  iP(t)  ^ 

dt  dt 

V  ^an  Cdn  j 

V^an  ^dn  y 

(44) 


where  kn  =  1  if  plasma  is  ignited  in  the  circuit  and  zero  otherwise,  and  the  plasma  current  is  given  by 

Ir(t)  =  ^L[vapp(t)-Vn(t)] 


(45) 


The  resistance  takes  on  values  Rnf  or  Rnb  depending  on  the  direction  of  the  current  in  the  plasma. 
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Equation  (44)  is  solved  using  multistage  Runge-Kutta  time  integration.  Then,  the  applied  voltage  and 
the  subcircuit  voltages  are  used  as  boundary  conditions  in  a  two-dimensional  solution  to  the  Poisson  equation, 
which  describes  the  electric  potential  4>: 


a2( |> 

0x2  +  0y2  ~  X2D 


(46) 


The  term  '/.[>  is  the  characteristic  length  for  electrostatic  shielding  in  plasma,  called  the  Debye  length. 

The  solution  domain  for  Eq  (46)  is  a  section  extracted  from  the  overall  Navier-Stokes  grid  in  the 
immediate  vicinity  of  the  plasma  actuator.  This  grid  section  is  subdivided  using  a  nonlinear  LaGrange 
interpolation  procedure12.  Once  the  body  forces  have  been  determined,  they  are  mapped  back  to  the  Navier- 
Stokes  grid.  This  procedure  allows  the  plasma  actuator  to  be  modeled  in  fine  detail  without  requiring  the  flow 
solver  to  operate  on  an  inordinate  number  of  field  points.  The  boundary  conditions  are  applied  to  the  domain 
boundaries  as  shown  in  Figure  4. 


Remainder  of  Boundary 
0  =  0 


Exposed  Insulated 
Electrode  Electrode 

o  =  vapp  o  =  vn(t) 


Figure  4:  Poisson  Equation  Solution  Domain  with  Boundary  Conditions 


In  practice,  the  Poisson  equation  is  applied  only  to  the  region  above  the  insulated  electrode  where 
plasma  is  present.  Elsewhere,  the  LaPlace  equation  is  solved: 


4+^i  =  0 

0x2  dy~ 


(47) 


Once  Eqs  (46)  and  (47)  are  solved,  the  body  force  per  volume  of  plasma  at  each  point  in  the  solution  domain 
can  be  computed.  The  electric  field  E  is  calculated  from  the  electric  potential: 

E  =  -V(|>  (48) 


And  the  body  force  is: 


fb=-l 


<t>E  = 


V^D  J 


{j>V<t> 


(49) 


Since  the  computational  grids  used  for  real  configurations  such  as  airfoils  are  not  Cartesian  with  even 
spacing  in  both  coordinate  directions,  the  Poisson/Laplace  equations  are  recast  from  Cartesian  coordinates  (x,y) 


12 


to  general  curvilinear  coordinates  (c,r|)  in  the  flow  solver.  The  curvilinear  coordinate  system  is  Cartesian  with 
uniform  spacing  of  one  unit  between  points  in  both  \  and  r|. 

The  transformed  Poisson  equation  is13 

J2 (a^  -  2b<\>^  +  c<|)ntl  +  dc^  +  e<j>,= )  =  -f-  (50) 

A,d 


where 


J=  1 

x§y,l  -YtjXr, 

(51) 

a  —  x“  +  y^ 
b  =  x^x^  +  y^yn 

9  9 

c  =  x^  +yZ 
d  =  j(cxy,  -Px?) 

(52) 

e  =  jfPx^  -ayj 

a  =  axEt  -  2bx=  +  ex 

(53) 

P  =  ayK-2by§tl+cyriT1 

(54) 

effect  the  coordinate  transformation. 

Eq  (50)  is  solved  using  a  finite  difference,  successive  line  over-relaxation  (SLOR)  algorithm.  Using 


the  standard  finite  difference  formulae9 


du  „  ui+i  -u,  i 


dx 


2  Ax 


O 


(Axf 


du  ui+1  —  2iij  +  2ui_1 


ax" 


(Ax)2 


-O 


(ax)2; 


(55) 


<3  U  ui+l,j+l  ui-l,j+l  ui+l,j-l  +  ui-l,. 


dxdy 


4AxAy 


^i  +  o[(Ax)2,(Ay)2 


the  discrete  form  of  Eq  (50)  for  solution  at  the  ith  grid  line  is 


'  ,._V 

C'j  2 

v  z  j 


TH 


]2X2 


r^-  +  2(aij  +cij) 


( 


C'j  +  2 

v  z  j 


h.j+1  _  f^ij 


(56) 


where 


(57) 


using  the  currently  available  values  of  <|).  Eq  (57)  is  solved  using  the  Thomas  algorithm  for  tridiagonal 
matrices6. 
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The  electric  field  components  [Eq  (48)]  are  computed  using  the  following  reverse  coordinate 
transformations: 


Ex  =  “<l>x  =-fex(l)^  +  11x<l>J 

Ey  =-<|>y  =  -fey<t>i;+Ry<t>J 


(58) 


where 


^x  =  Jy^  =  -Jx 

Tlx  =  -Jy^  fly  = 


(59) 


In  the  plasma  actuator  model  supplied  by  the  University  of  Notre  Dame,  the  sign  of  the  body  force  is 
the  opposite  of  that  shown  in  Eq  (49).  The  reason  for  the  body  force  sign  reversal  in  the  Notre  Dame  code  is 
unknown;  this  issue  had  not  been  resolved  with  the  Notre  Dame  researchers  at  the  time  of  this  writing. 


RESULTS 

Two  example  calculations  are  presented  here.  The  first  is  an  airfoil  flowfield  solution,  calculated 
without  the  plasma  actuator  model,  which  demonstrates  the  accuracy  of  the  flow  solver  and  illustrates  the 
importance  of  modeling  turbulence  in  the  flow  solutions.  The  second  solution  employs  a  plasma  actuator  to 
induce  motion  in  quiescent  air  over  a  flat  plate. 

RAE  2822  Airfoil 

Laminar  and  turbulent  solutions  for  a  standard  RAE  2822  airfoil  CFD  code  validation  case  (Mach 
0.729,  a  =  2.31°)  are  contrasted  in  Figure  5.  The  surface  pressure  distributions  from  the  laminar  and  turbulent 
solutions  are  compared  to  wind  tunnel  test  data  for  this  case  in  Figure  6.  Agreement  of  the  turbulent  solution 
pressure  with  the  data  is  excellent. 

These  results  illustrate  the  importance  of  turbulence  modeling  for  accurate  flow  characterization. 
Plasma  actuators  are  being  investigated  for  use  in  fluid  separation  control.  It  is  imperative  that  tools  developed 
to  simulate  plasma  actuator  effects  adequately  model  the  basic  fluid  flowfield  so  that  valid  conclusions  can  be 
drawn  from  the  simulations  about  the  fluid  flow  manipulation  capability  of  these  devices.  Separation  should 
only  be  indicated  in  the  flow  simulations  in  those  cases  in  which  separation  is  likely  to  occur.  In  the  present 
case,  the  laminar  solution  inaccurately  produced  large  scale  flow  separation  and  grossly  inaccurate  pressure 
fields.  Including  turbulence  modeling  in  the  simulation  eliminated  the  separation  and  greatly  improved  the 
agreement  between  the  simulation  and  test  data. 
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(b)  Turbulent 

Figure  5:  Laminar  and  Turbulent  Solutions  at  Corresponding  Solution  Times  for  the  RAE  2822  Airfoil, 

Mach  0.729,  a  =  2.31° 


Figure  6:  RAE  2822  Airfoil  Pressure  Data  (Points),  Laminar  Simulation  (Red  Line)  and  Turbulent 

Simulation  (Blue  Line) 
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Flat  Plate  with  Plasma  Actuator 


The  plasma  actuator  model  was  tested  by  evaluating  the  effect  of  a  plasma  actuator  mounted  on  a  flat 
plate  on  the  surrounding  quiescent  air.  The  solution  domain  for  this  case  was  very  similar  to  the  domain 
illustrated  in  Figure  4.  Both  the  exposed  and  hidden  electrodes  had  a  length  of  Vi  inch  (0.0127  m).  The 
dielectric 

material  was  3-mil  thick  (7.62x  10"5  m)  Kapton  film  with  a  dielectric  coefficient  of  3.5.  The  voltage  of  5  kV  was 
applied  at  an  AC  frequency  of  2  kHz. 

The  2m/  1  m.  161x81  point  computational  grid  was  clustered  near  the  plasma  actuator,  which  was 
centered  on  the  lower  boundary  of  the  solution  domain.  There  were  21  evenly  spaced  points  on  each  electrode. 
The  spacing  normal  to  the  wall  was  a  very  coarse  10"4  m.  The  computational  grid  is  shown  in  Figure  7. 


(a)  Solution  Domain  (b)  Plasma  Actuator  Region 

Figure  7:  Computational  Grid  for  Plasma  Actuator  on  a  Flat  Plate 


The  flat  plate  solution  was  ran  using  a  computational  time  step  of  2  ps  for  4,000  time  steps.  The 
development  of  a  left-to-right  near- wall  flow  is  shown  in  Figure  8  for  the  first  two  AC  cycles.  Note  the 
development  and  convection  of  the  large  initial  starting  vortex,  and  the  smaller  secondary  vortex  that  forms 
during  the  second  AC  cycle.  The  extent  of  induced  fluid  motion  at  different  simulation  times  is  shown  in  Figure 
9.  The  simulated  fluid  motion  is  consistent  with  the  experimental  observations  of  Enloe  et  al'. 

CONCLUSION  AND  RECOMMENDATIONS 

A  time-accurate  Navier-Stokes  flow  solver  with  an  integrated  plasma  actuator  body  force  model  has  been 
developed.  The  base  flow  solver’s  accuracy  has  been  demonstrated  using  standard  validation  cases.  Operation 
of  the  plasma  actuator  model  has  also  been  successfully  demonstrated,  although  a  thorough  validation  study  has 
not  been  conducted. 

Since  AFOSR  is  no  longer  providing  funding  to  AFRL/MNAC  for  this  research,  AFRL/MNAC’s 
research  effort  is  being  discontinued.  However,  under  an  agreement  reached  with  the  University  of  Notre  Dame, 
the  code  will  be  provided  to  researchers  there  in  the  hope  that  it  can  be  fully  developed  into  a  useful  analysis 
tool.  The  following  suggestions  are  offered  for  future  development  of  the  code: 
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(i)  t  =  1.8  ms  (j)t  =  2.0ms 

Figure  8:  Near-Wall  Momentum  Vectors  During  the  First  Two  Plasma  Actuator  AC  Cycles 


l»wv" 


(d)  t  =  8  ms 

Figure  9:  X-Momentum  Propagation  During  the  Plasma  Actuator  Simulation 

❖  Thorough  Validation  of  the  Plasma  Actuator  Model.  Although  great  care  was  taken  in  the  plasma 
actuator  model  implementation,  the  model  originators  should  carefully  examine  the  code  and  output  to 
ensure  the  proper  implementation  and  operation  of  the  plasma  actuator  model. 

❖  Multiple  Boundary  Conditions  on  Grid  Boundaries:  The  flow  solver  currently  allows  the  use  of  only  one 
boundary  condition  per  grid  boundary.  This  limits  the  allowable  computational  grid  types  to  H-  and  O- 
topologies.  This  simple  change  would  allow  the  use  of  C-topology  grids,  which  are  better  suited  for  airfoil 
analysis  than  the  presently  allowed  O-topology  grids.  This  capability  was  included  in  a  three-dimensional 
flow  solver  developed  for  this  research  effort  and  should  serve  as  an  example  for  implementation  of  a 
similar  capability  in  the  two-dimensional  code. 

❖  Characteristic  Outer  Boundary  Condition:  Currently,  only  extrapolation  and  freestream  conditions  can  be 
imposed  at  the  outer  boundary.  This  introduces  a  requirement  that  the  outer  grid  boundary  of  an  airfoil  grid 
must  be  a  great  physical  distance  from  the  inner  grid  boundary  in  order  to  adequately  simulate  real 
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flowfields.  This  capability  was  also  included  in  the  previously  mentioned  three-dimensional  flow  solver 
and  can  be  easily  implemented  in  the  two-dimensional  code. 

Multiple  Plasma  Actuators'.  Due  to  time  constraints,  this  capability  was  not  written  into  the  code. 
However,  such  a  capability  is  essential  if  this  analysis  capability  is  to  be  maximized. 

Parallel  Processing  Capability.  While  the  code  runs  rapidly  on  a  single  CPU,  analyses  could  be  conducted 
much  more  rapidly  if  this  capability  was  added. 
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